stanEx5
fn=function(mdl,data,smp=500,wrm=100){
mle=mdl$optimize(data=data)
print(mle)
mcmc=goMCMC(mdl,data,smp,wrm)
mcmc$metadata()$stan_variables
seeMCMC(mcmc,'m')
y0=mcmc$draws('y0')
smy0=summary(y0)
grid.arrange(
qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))
par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')
m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)
xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)
qplot(x,y,col=I('red'))+
geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=m),xy,col='black')
}ex12
hierarchical model
class c=1~k
b0c~N(b00,sb0)
b1c~N(b10,sb1)
yc~N(b0c+b1c*x,s)n=50
k=9
x=runif(n,0,20)
c=sample(letters[1:k],n,replace=T)
b00=rnorm(k,10,5)
b0=b00[factor(c)]
b10=rnorm(k,2,1)
b1=b10[factor(c)]
y=rnorm(n,b0+b1*x,5)
qplot(x,y,shape=c,size=I(2))+
scale_shape_manual(values=1:k)
qplot(x,y,col=c)
n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcitionestimate as no class
y~N(b00+b10*x,s)ex8-0.stan
data{
int N;
vector[N] x;
vector[N] y;
int N1;
vector[N1] x1;
}
parameters{
real b0;
real b1;
real<lower=0> s;
}
model{
y~normal(b0+b1*x,s);
}
generated quantities{
vector[N] m0;
vector[N] y0;
for(i in 1:N){
m0[i]=b0+b1*x[i];
y0[i]=normal_rng(m0[i],s);
}
vector[N1] m1;
vector[N1] y1;
for(i in 1:N1){
m1[i]=b0+b1*x1[i];
y1[i]=normal_rng(m1[i],s);
}
}mdl=cmdstan_model('./ex8-0.stan')
data=list(N=n,x=x,y=y,N1=n1,x1=x1)
fn(mdl,data)## Initial log joint probability = -7288.33
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 20 -135.446 0.000294158 0.0039898 1 1 35
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.2 seconds.
## variable estimate
## lp__ -135.45
## b0 7.07
## b1 1.92
## s 9.11
## m0[1] 14.63
## m0[2] 15.39
## m0[3] 39.30
## m0[4] 16.67
## m0[5] 44.47
## m0[6] 44.65
## m0[7] 29.04
## m0[8] 13.85
## m0[9] 38.77
## m0[10] 29.90
## m0[11] 8.42
## m0[12] 44.15
## m0[13] 44.32
## m0[14] 33.87
## m0[15] 24.36
## m0[16] 26.71
## m0[17] 17.59
## m0[18] 20.37
## m0[19] 10.98
## m0[20] 15.36
## m0[21] 11.53
## m0[22] 36.15
## m0[23] 40.70
## m0[24] 15.24
## m0[25] 23.22
## m0[26] 18.09
## m0[27] 36.32
## m0[28] 23.05
## m0[29] 30.48
## m0[30] 36.81
## m0[31] 19.11
## m0[32] 22.81
## m0[33] 17.16
## m0[34] 43.25
## m0[35] 35.80
## m0[36] 19.27
## m0[37] 19.67
## m0[38] 37.20
## m0[39] 32.30
## m0[40] 34.85
## m0[41] 31.38
## m0[42] 41.09
## m0[43] 13.92
## m0[44] 27.25
## m0[45] 40.50
## m0[46] 41.42
## m0[47] 25.03
## m0[48] 20.34
## m0[49] 19.14
## m0[50] 41.19
## y0[1] 12.35
## y0[2] 7.84
## y0[3] 40.86
## y0[4] 19.76
## y0[5] 57.44
## y0[6] 48.30
## y0[7] 39.57
## y0[8] 14.43
## y0[9] 30.93
## y0[10] 28.43
## y0[11] -5.66
## y0[12] 41.27
## y0[13] 43.45
## y0[14] 34.74
## y0[15] 23.13
## y0[16] 9.36
## y0[17] 17.71
## y0[18] 30.73
## y0[19] 5.88
## y0[20] 7.10
## y0[21] -0.90
## y0[22] 38.26
## y0[23] 44.15
## y0[24] 3.58
## y0[25] 13.92
## y0[26] 26.23
## y0[27] 38.65
## y0[28] 10.51
## y0[29] 22.36
## y0[30] 34.38
## y0[31] 45.13
## y0[32] 28.55
## y0[33] 8.98
## y0[34] 54.68
## y0[35] 39.62
## y0[36] 27.56
## y0[37] 18.59
## y0[38] 31.25
## y0[39] 24.49
## y0[40] 37.17
## y0[41] 35.15
## y0[42] 47.67
## y0[43] 23.43
## y0[44] 25.27
## y0[45] 25.17
## y0[46] 46.84
## y0[47] 17.49
## y0[48] 31.01
## y0[49] 22.22
## y0[50] 36.26
## m1[1] 8.42
## m1[2] 12.45
## m1[3] 16.47
## m1[4] 20.50
## m1[5] 24.52
## m1[6] 28.55
## m1[7] 32.57
## m1[8] 36.60
## m1[9] 40.62
## m1[10] 44.65
## y1[1] 1.69
## y1[2] -3.07
## y1[3] 13.57
## y1[4] 17.58
## y1[5] 27.78
## y1[6] 31.94
## y1[7] 20.63
## y1[8] 48.37
## y1[9] 26.64
## y1[10] 50.46
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.3 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -134.70 -134.43 1.18 0.96 -136.92 -133.39 1.00 654 1121
## b0 7.27 7.27 2.65 2.73 2.96 11.49 1.01 353 302
## b1 1.90 1.90 0.22 0.23 1.52 2.25 1.00 435 606
## s 9.54 9.43 1.02 0.99 8.04 11.30 1.00 1851 1491
## m0[1] 14.76 14.77 1.94 2.01 11.61 17.93 1.01 372 339
## m0[2] 15.51 15.56 1.88 1.94 12.45 18.58 1.01 377 344
## m0[3] 39.19 39.22 1.98 1.99 35.88 42.33 1.00 1414 1249
## m0[4] 16.77 16.80 1.77 1.81 13.86 19.69 1.01 389 344
## m0[5] 44.31 44.31 2.46 2.47 40.30 48.34 1.00 994 1094
## m0[6] 44.49 44.49 2.48 2.50 40.46 48.56 1.00 984 1106
## m0[7] 29.04 29.06 1.35 1.34 26.79 31.24 1.00 1560 1204
## m0[8] 13.99 14.00 2.01 2.07 10.74 17.23 1.01 367 337
## m0[9] 38.67 38.70 1.93 1.94 35.48 41.74 1.00 1474 1259
## m0[10] 29.88 29.90 1.38 1.37 27.60 32.10 1.00 1763 1238
## m0[11] 8.61 8.60 2.52 2.57 4.55 12.62 1.01 354 296
## m0[12] 44.00 44.01 2.43 2.44 40.03 47.98 1.00 1011 1083
## m0[13] 44.17 44.17 2.45 2.46 40.18 48.18 1.00 1001 1094
## m0[14] 33.81 33.85 1.57 1.61 31.25 36.40 1.00 2017 1399
## m0[15] 24.40 24.39 1.36 1.31 22.18 26.68 1.00 684 815
## m0[16] 26.72 26.73 1.33 1.30 24.53 28.94 1.00 1011 961
## m0[17] 17.69 17.71 1.70 1.73 14.90 20.48 1.01 400 358
## m0[18] 20.44 20.45 1.52 1.51 17.87 22.97 1.00 457 447
## m0[19] 11.14 11.15 2.27 2.31 7.49 14.75 1.01 357 311
## m0[20] 15.48 15.53 1.88 1.94 12.42 18.55 1.01 377 349
## m0[21] 11.68 11.69 2.22 2.27 8.12 15.21 1.01 358 305
## m0[22] 36.07 36.13 1.72 1.76 33.25 38.83 1.00 1752 1376
## m0[23] 40.59 40.59 2.11 2.12 37.11 43.97 1.00 1260 1177
## m0[24] 15.36 15.40 1.89 1.95 12.30 18.45 1.01 376 344
## m0[25] 23.27 23.28 1.39 1.36 20.95 25.64 1.00 589 720
## m0[26] 18.19 18.20 1.67 1.67 15.43 20.95 1.01 407 396
## m0[27] 36.25 36.31 1.74 1.77 33.41 39.01 1.00 1734 1276
## m0[28] 23.10 23.11 1.40 1.38 20.76 25.48 1.00 578 719
## m0[29] 30.45 30.46 1.40 1.41 28.15 32.72 1.00 1793 1254
## m0[30] 36.73 36.80 1.78 1.79 33.84 39.56 1.00 1677 1319
## m0[31] 19.20 19.19 1.60 1.60 16.54 21.80 1.01 426 412
## m0[32] 22.86 22.87 1.40 1.39 20.50 25.26 1.00 563 703
## m0[33] 17.27 17.28 1.74 1.78 14.43 20.11 1.01 395 358
## m0[34] 43.10 43.11 2.34 2.37 39.24 46.94 1.00 1065 1082
## m0[35] 35.73 35.78 1.70 1.74 32.95 38.45 1.00 1791 1354
## m0[36] 19.35 19.35 1.59 1.58 16.71 21.95 1.01 429 423
## m0[37] 19.75 19.76 1.56 1.56 17.13 22.34 1.01 439 432
## m0[38] 37.12 37.17 1.81 1.81 34.16 40.00 1.00 1678 1302
## m0[39] 32.26 32.28 1.48 1.52 29.80 34.67 1.00 2022 1413
## m0[40] 34.79 34.83 1.63 1.69 32.14 37.42 1.00 1888 1422
## m0[41] 31.35 31.35 1.43 1.46 28.99 33.69 1.00 1910 1288
## m0[42] 40.97 40.98 2.14 2.15 37.45 44.41 1.00 1226 1177
## m0[43] 14.05 14.07 2.00 2.07 10.82 17.29 1.01 367 337
## m0[44] 27.25 27.26 1.33 1.30 25.07 29.46 1.00 1121 1083
## m0[45] 40.39 40.40 2.09 2.11 36.94 43.74 1.00 1294 1177
## m0[46] 41.29 41.30 2.17 2.19 37.72 44.78 1.00 1199 1177
## m0[47] 25.06 25.05 1.34 1.31 22.85 27.33 1.00 755 875
## m0[48] 20.42 20.43 1.52 1.51 17.84 22.95 1.00 456 447
## m0[49] 19.23 19.23 1.60 1.59 16.58 21.83 1.01 426 412
## m0[50] 41.07 41.08 2.15 2.16 37.53 44.52 1.00 1217 1177
## y0[1] 14.68 15.02 9.79 9.53 -1.73 30.53 1.00 1621 1929
## y0[2] 15.05 15.16 9.65 9.51 -0.69 30.65 1.00 1965 1745
## y0[3] 39.40 39.43 9.59 9.26 23.70 55.27 1.00 1975 1636
## y0[4] 16.81 16.82 9.87 9.70 0.76 32.58 1.00 1823 1853
## y0[5] 44.67 45.01 10.13 10.08 27.97 60.55 1.00 1899 2013
## y0[6] 44.64 44.56 9.88 9.68 28.74 60.54 1.00 1778 1625
## y0[7] 29.20 29.06 9.56 9.33 13.46 45.18 1.00 2004 1993
## y0[8] 14.53 14.83 9.81 9.82 -1.54 31.22 1.00 1718 1730
## y0[9] 38.93 39.14 9.80 9.91 23.28 55.07 1.00 2074 2147
## y0[10] 30.01 30.33 9.54 9.16 14.32 45.85 1.00 2121 1897
## y0[11] 8.88 8.83 10.01 10.16 -7.32 25.09 1.00 1425 1765
## y0[12] 44.11 43.86 9.78 9.63 28.26 60.57 1.00 1703 1968
## y0[13] 44.60 44.48 10.03 10.00 28.72 61.03 1.00 1907 1914
## y0[14] 33.62 33.67 9.46 9.52 18.11 48.91 1.00 2084 1776
## y0[15] 24.32 24.40 9.75 9.29 8.05 40.21 1.00 2011 1927
## y0[16] 26.52 26.66 9.89 9.71 10.17 42.52 1.00 1985 1895
## y0[17] 17.40 17.49 9.40 9.11 2.24 33.05 1.00 1467 1856
## y0[18] 20.17 20.37 9.76 9.50 3.65 35.48 1.00 1773 1918
## y0[19] 10.92 11.05 10.03 9.74 -6.40 27.30 1.00 1605 1822
## y0[20] 15.22 15.09 9.62 9.28 -0.21 31.71 1.00 2012 2040
## y0[21] 11.84 11.75 9.78 9.40 -4.31 27.91 1.00 2007 2022
## y0[22] 35.99 36.04 9.85 9.54 19.75 52.74 1.00 1839 1966
## y0[23] 40.30 40.20 9.73 9.50 24.32 56.55 1.00 1850 2101
## y0[24] 14.98 15.16 9.78 9.72 -1.03 30.74 1.00 1740 1933
## y0[25] 23.10 23.08 9.69 9.65 7.18 38.83 1.00 1985 1894
## y0[26] 17.90 17.87 9.56 9.44 2.39 33.49 1.00 1725 1820
## y0[27] 36.20 35.89 10.09 9.99 20.17 52.61 1.00 1887 1857
## y0[28] 23.20 23.17 9.79 9.75 7.14 39.26 1.00 2120 1862
## y0[29] 30.48 30.58 9.69 9.59 14.40 46.51 1.00 1978 2057
## y0[30] 36.49 36.71 9.74 9.52 19.77 52.15 1.00 1882 1834
## y0[31] 19.30 19.20 9.85 9.73 3.22 35.19 1.00 1874 2056
## y0[32] 22.49 22.65 9.68 9.29 6.64 38.45 1.00 1540 1815
## y0[33] 17.16 17.12 9.71 9.74 1.61 33.53 1.00 1487 1877
## y0[34] 43.20 42.82 9.93 9.95 27.64 59.43 1.00 1841 2053
## y0[35] 35.64 36.00 9.82 10.04 19.41 51.69 1.00 1877 1928
## y0[36] 18.87 18.75 9.46 9.19 3.48 34.64 1.00 1884 2007
## y0[37] 19.68 20.00 9.65 9.67 3.59 35.53 1.00 2008 1930
## y0[38] 37.22 37.16 9.75 9.58 21.44 53.11 1.00 2254 1746
## y0[39] 32.45 32.55 9.47 9.06 16.60 47.96 1.00 2158 2028
## y0[40] 34.78 34.59 10.04 9.82 18.71 51.53 1.00 2179 1913
## y0[41] 31.37 31.31 9.79 9.87 15.72 47.26 1.00 2046 1826
## y0[42] 41.12 41.15 9.96 9.85 24.75 57.09 1.00 1966 1957
## y0[43] 13.91 14.12 9.75 9.50 -1.81 29.72 1.00 1893 1929
## y0[44] 27.41 27.43 9.72 9.68 11.01 43.45 1.00 2043 1957
## y0[45] 40.08 39.67 9.84 9.63 24.56 56.34 1.00 1896 2101
## y0[46] 40.98 40.70 9.51 9.17 25.41 57.13 1.01 1803 2056
## y0[47] 25.19 25.18 9.68 9.65 9.21 40.73 1.00 1930 1845
## y0[48] 20.12 20.13 9.69 9.52 4.25 36.20 1.00 1914 2003
## y0[49] 19.14 19.28 9.71 9.85 2.66 35.18 1.00 1848 2053
## y0[50] 41.18 41.06 9.84 9.92 25.22 57.56 1.00 1907 1636
## m1[1] 8.61 8.60 2.52 2.57 4.55 12.62 1.01 354 296
## m1[2] 12.60 12.62 2.13 2.19 9.15 16.04 1.01 361 310
## m1[3] 16.58 16.62 1.79 1.82 13.64 19.53 1.01 387 344
## m1[4] 20.57 20.58 1.51 1.50 18.01 23.09 1.00 461 450
## m1[5] 24.56 24.55 1.35 1.32 22.34 26.83 1.00 700 802
## m1[6] 28.55 28.56 1.34 1.33 26.32 30.75 1.00 1438 1209
## m1[7] 32.53 32.53 1.49 1.53 30.05 34.96 1.00 2025 1422
## m1[8] 36.52 36.60 1.76 1.79 33.65 39.32 1.00 1701 1276
## m1[9] 40.51 40.51 2.10 2.12 37.04 43.88 1.00 1283 1177
## m1[10] 44.49 44.49 2.48 2.50 40.46 48.56 1.00 984 1106
## y1[1] 8.90 8.97 10.10 9.99 -7.86 25.42 1.00 1445 1904
## y1[2] 12.67 12.83 9.86 9.54 -3.12 28.87 1.00 1660 1618
## y1[3] 16.78 16.90 9.85 9.65 0.21 32.38 1.00 1888 1563
## y1[4] 20.53 20.20 9.79 9.65 4.45 36.95 1.00 1799 1936
## y1[5] 24.22 24.12 9.68 9.40 8.24 40.32 1.00 1875 1854
## y1[6] 28.47 28.22 9.73 9.96 12.88 44.60 1.00 2045 1845
## y1[7] 32.72 32.69 9.54 9.57 16.36 48.27 1.00 1859 1974
## y1[8] 36.74 36.83 9.57 9.49 21.17 52.02 1.00 1462 2097
## y1[9] 40.22 40.23 10.00 10.23 23.96 56.12 1.00 1823 1750
## y1[10] 44.40 44.12 9.71 9.45 28.55 60.17 1.00 1662 1933
estimate as independent class
all b0l,b1l do not have a distribution
class c=1~k
yc~N(b0c+b1c*x,s) ex12-1.stan
data{
int N;
int K;
vector[N] x;
vector[N] y;
array[N] int c;
}
parameters{
vector[K] b0;
vector[K] b1;
real<lower=0> s;
}
model{
for(i in 1:N)
y[i]~normal(b0[c[i]]+b1[c[i]]*x[i],s);
}
generated quantities{
vector[N] m0;
vector[N] y0;
for(i in 1:N){
m0[i]=b0[c[i]]+b1[c[i]]*x[i];
y0[i]=normal_rng(m0[i],s);
}
}mdl=cmdstan_model('./ex12-1.stan')
data=list(N=n,K=k,x=x,y=y,c=factor(c))
mle=mdl$optimize(data=data)## Initial log joint probability = -1.23086e+06
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## Error evaluating model log probability: Non-finite gradient.
## Error evaluating model log probability: Non-finite gradient.
## Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/tmp/Rtmpwb3xsj/model-20a93d7d50fa.stan', line 15, column 4 to column 42)
## Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/tmp/Rtmpwb3xsj/model-20a93d7d50fa.stan', line 15, column 4 to column 42)
## Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/tmp/Rtmpwb3xsj/model-20a93d7d50fa.stan', line 15, column 4 to column 42)
## Error evaluating model log probability: Non-finite gradient.
## 99 -87.0685 0.426789 1.49779 1 1 166
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 199 -86.9818 0.00250448 0.0405116 0.6326 0.6326 275
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 222 -86.9817 0.000339057 0.00138477 1 1 300
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -86.98
## b0[1] -1.59
## b0[2] 15.59
## b0[3] 21.29
## b0[4] 5.07
## b0[5] 5.36
## b0[6] 3.55
## b0[7] 8.14
## b0[8] -2.70
## b0[9] -1.45
## b1[1] 3.85
## b1[2] 1.03
## b1[3] 1.70
## b1[4] 0.18
## b1[5] 2.14
## b1[6] 2.27
## b1[7] 1.64
## b1[8] 1.55
## b1[9] 3.17
## s 3.45
## m0[1] 12.49
## m0[2] 14.65
## m0[3] 32.99
## m0[4] 29.80
## m0[5] 47.14
## m0[6] 35.88
## m0[7] 15.03
## m0[8] 27.30
## m0[9] 49.41
## m0[10] 44.29
## m0[11] 5.15
## m0[12] 47.43
## m0[13] 46.98
## m0[14] 18.92
## m0[15] 24.01
## m0[16] 27.30
## m0[17] 15.98
## m0[18] 8.03
## m0[19] 8.18
## m0[20] 20.06
## m0[21] 5.48
## m0[22] 20.76
## m0[23] 24.44
## m0[24] 15.15
## m0[25] 35.61
## m0[26] 16.59
## m0[27] 47.24
## m0[28] 22.46
## m0[29] 37.34
## m0[30] 38.59
## m0[31] 18.50
## m0[32] 35.25
## m0[33] 5.44
## m0[34] 53.38
## m0[35] 31.10
## m0[36] 6.19
## m0[37] 22.39
## m0[38] 39.03
## m0[39] 33.55
## m0[40] 31.98
## m0[41] 32.52
## m0[42] 51.47
## m0[43] 11.65
## m0[44] 27.42
## m0[45] 33.64
## m0[46] 37.61
## m0[47] 37.22
## m0[48] 33.06
## m0[49] 7.04
## m0[50] 43.92
## y0[1] 15.49
## y0[2] 15.04
## y0[3] 36.53
## y0[4] 36.96
## y0[5] 49.63
## y0[6] 33.23
## y0[7] 14.31
## y0[8] 24.11
## y0[9] 48.09
## y0[10] 46.05
## y0[11] 4.25
## y0[12] 50.60
## y0[13] 46.70
## y0[14] 19.58
## y0[15] 22.96
## y0[16] 30.07
## y0[17] 16.09
## y0[18] 8.47
## y0[19] 11.16
## y0[20] 16.61
## y0[21] 4.92
## y0[22] 26.14
## y0[23] 30.23
## y0[24] 10.91
## y0[25] 38.87
## y0[26] 12.09
## y0[27] 47.69
## y0[28] 20.29
## y0[29] 44.66
## y0[30] 41.58
## y0[31] 16.53
## y0[32] 36.94
## y0[33] 4.66
## y0[34] 57.22
## y0[35] 27.04
## y0[36] 5.57
## y0[37] 19.46
## y0[38] 39.37
## y0[39] 38.38
## y0[40] 31.87
## y0[41] 32.49
## y0[42] 51.40
## y0[43] 16.02
## y0[44] 28.10
## y0[45] 30.18
## y0[46] 45.54
## y0[47] 35.36
## y0[48] 40.43
## y0[49] 3.77
## y0[50] 44.06mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 2.4 seconds.
## Chain 2 finished in 2.4 seconds.
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 finished in 2.5 seconds.
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 finished in 2.6 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 2.5 seconds.
## Total execution time: 2.7 seconds.seeMCMC(mcmc,'m')## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -97.13 -96.75 3.75 3.65 -103.98 -91.71 1.01 653 1058
## b0[1] -333.81 0.10 729.64 372.80 -1710.25 379.46 1.88 5 28
## b0[2] 15.71 15.82 4.52 4.37 8.11 23.35 1.00 2422 1268
## b0[3] 21.20 21.26 3.15 3.08 16.07 26.12 1.00 3419 1349
## b0[4] 5.07 4.87 7.59 7.32 -7.04 17.87 1.00 1171 1040
## b0[5] 5.42 5.59 4.90 4.84 -2.85 13.28 1.01 2325 1464
## b0[6] 3.58 3.59 2.37 2.33 -0.11 7.53 1.01 2681 923
## b0[7] 8.17 8.08 6.32 6.28 -2.23 18.96 1.01 1957 1206
## b0[8] -2.66 -2.68 4.49 4.40 -9.93 4.68 1.00 2217 1172
## b0[9] -1.40 -1.74 7.34 7.42 -13.18 10.80 1.00 1379 1164
## b1[1] 31.73 3.79 61.22 31.14 -28.31 147.29 1.88 5 27
## b1[2] 1.03 1.03 0.31 0.30 0.50 1.55 1.00 2208 1466
## b1[3] 1.71 1.70 0.26 0.27 1.29 2.12 1.00 3062 1560
## b1[4] 0.14 0.18 1.57 1.48 -2.53 2.56 1.01 1208 1130
## b1[5] 2.14 2.13 0.34 0.33 1.58 2.70 1.00 2358 1338
## b1[6] 2.27 2.27 0.23 0.23 1.89 2.63 1.00 2213 1566
## b1[7] 1.63 1.64 0.47 0.46 0.88 2.39 1.01 1750 1379
## b1[8] 1.54 1.54 0.38 0.37 0.92 2.17 1.00 2291 1259
## b1[9] 3.17 3.18 0.85 0.85 1.76 4.57 1.00 1548 1251
## s 4.44 4.38 0.57 0.54 3.62 5.47 1.00 1200 1366
## m0[1] 12.52 12.53 1.74 1.70 9.67 15.34 1.00 2559 1198
## m0[2] 14.69 14.82 3.55 3.51 8.74 20.37 1.00 2314 1339
## m0[3] 33.00 33.00 2.15 2.00 29.53 36.55 1.01 1830 989
## m0[4] 29.75 29.82 2.08 2.04 26.38 33.04 1.00 3134 1408
## m0[5] 47.12 47.11 2.53 2.45 43.10 51.29 1.00 2200 1276
## m0[6] 35.87 35.87 2.71 2.61 31.45 40.30 1.01 1840 1047
## m0[7] 15.03 15.06 1.69 1.63 12.25 17.84 1.00 2363 1303
## m0[8] 27.24 27.31 2.37 2.34 23.39 30.93 1.00 3277 1371
## m0[9] 49.43 49.40 2.07 2.04 46.01 52.82 1.00 2173 1445
## m0[10] 44.41 44.52 4.43 4.06 37.05 51.77 1.00 2402 1261
## m0[11] 5.18 5.19 2.24 2.17 1.65 8.86 1.01 2700 942
## m0[12] 47.45 47.44 3.02 2.98 42.57 52.33 1.00 1986 1567
## m0[13] 46.96 46.95 2.51 2.43 42.98 51.11 1.00 2198 1276
## m0[14] 18.91 18.92 2.05 2.00 15.45 22.26 1.00 2295 1404
## m0[15] 24.03 24.02 1.49 1.47 21.61 26.42 1.00 2085 1807
## m0[16] 27.32 27.34 1.99 1.96 24.07 30.61 1.00 2224 1675
## m0[17] 16.03 15.98 3.40 3.52 10.53 21.72 1.00 1359 1331
## m0[18] 8.04 7.99 2.25 2.20 4.47 11.84 1.00 2303 1460
## m0[19] 8.21 8.21 2.02 1.96 5.01 11.40 1.00 2716 976
## m0[20] 20.16 20.28 3.33 3.24 14.64 25.76 1.00 2438 1226
## m0[21] 5.40 5.40 4.57 4.53 -1.99 13.23 1.00 1246 1275
## m0[22] 20.75 20.78 2.35 2.30 16.78 24.59 1.00 2296 1266
## m0[23] 24.42 24.44 3.06 2.95 19.30 29.34 1.00 2327 1141
## m0[24] 15.14 15.11 4.58 4.55 7.70 22.96 1.00 1990 1351
## m0[25] 35.59 35.64 1.57 1.49 33.01 38.13 1.00 2578 1589
## m0[26] 16.62 16.63 1.56 1.50 14.05 19.08 1.00 2337 1647
## m0[27] 47.26 47.24 1.85 1.80 44.19 50.28 1.00 2115 1419
## m0[28] 22.49 22.48 1.47 1.42 20.07 24.86 1.00 2144 1874
## m0[29] 37.40 37.42 4.39 4.30 30.11 44.56 1.00 1883 1386
## m0[30] 38.58 38.54 1.72 1.60 35.80 41.40 1.00 2142 1226
## m0[31] 18.56 18.51 3.01 3.12 13.54 23.48 1.00 1418 1346
## m0[32] 35.22 35.27 1.59 1.51 32.61 37.80 1.00 2612 1598
## m0[33] 5.46 5.38 2.72 2.66 1.00 9.92 1.00 2257 1370
## m0[34] 53.42 53.42 2.54 2.51 49.23 57.55 1.00 2292 1555
## m0[35] 31.13 31.12 1.92 1.82 28.05 34.29 1.01 1937 962
## m0[36] 5.97 6.12 4.57 4.51 -1.58 13.37 1.01 1444 1220
## m0[37] 22.47 22.54 2.77 2.69 17.99 27.08 1.00 2497 1237
## m0[38] 39.02 38.98 1.75 1.64 36.20 41.85 1.00 2140 1211
## m0[39] 33.55 33.55 1.62 1.54 30.88 36.25 1.00 2219 1535
## m0[40] 31.88 31.91 2.78 2.78 27.42 36.30 1.01 1935 1388
## m0[41] 32.53 32.54 1.65 1.56 29.85 35.30 1.00 2210 1587
## m0[42] 51.50 51.49 2.30 2.28 47.74 55.27 1.00 2239 1525
## m0[43] 11.68 11.70 1.79 1.75 8.78 14.56 1.00 2622 1164
## m0[44] 27.45 27.41 1.58 1.64 24.91 30.00 1.00 2042 1560
## m0[45] 33.65 33.69 2.26 2.13 30.02 37.40 1.01 1822 959
## m0[46] 37.48 37.47 3.66 3.61 31.28 43.25 1.01 1749 1412
## m0[47] 37.20 37.24 1.49 1.42 34.75 39.65 1.00 2429 1517
## m0[48] 33.02 33.07 1.76 1.67 30.15 35.86 1.00 2851 1351
## m0[49] 7.06 6.97 2.42 2.35 3.15 11.09 1.00 2292 1379
## m0[50] 43.95 43.94 2.71 2.67 39.60 48.31 1.00 1980 1545
## y0[1] 12.58 12.59 4.92 4.77 4.42 20.43 1.00 2274 2059
## y0[2] 14.79 14.78 5.68 5.56 5.09 23.90 1.00 2165 1926
## y0[3] 32.98 33.05 4.99 4.79 24.71 41.03 1.00 1785 1743
## y0[4] 29.74 29.71 4.80 4.53 21.71 37.69 1.00 2109 1881
## y0[5] 47.00 46.98 5.04 5.19 38.68 55.10 1.00 2102 1840
## y0[6] 35.89 35.78 5.26 5.15 27.38 44.81 1.00 1919 1810
## y0[7] 15.07 14.95 4.70 4.65 7.39 22.91 1.00 1882 1970
## y0[8] 27.27 27.38 4.98 4.81 19.17 35.38 1.00 1966 1826
## y0[9] 49.41 49.32 5.00 4.89 40.94 57.59 1.00 1905 2011
## y0[10] 44.33 44.47 6.34 6.04 33.97 54.55 1.00 2078 1732
## y0[11] 5.01 5.05 5.02 4.84 -3.22 13.50 1.00 1903 1669
## y0[12] 47.33 47.31 5.32 5.28 38.76 55.98 1.00 2159 1931
## y0[13] 47.02 46.90 5.09 5.08 38.52 55.36 1.00 2087 2039
## y0[14] 18.82 19.05 4.99 5.03 10.71 26.69 1.00 1812 1883
## y0[15] 23.83 23.99 4.65 4.46 16.06 31.45 1.00 1620 1820
## y0[16] 27.36 27.28 4.89 4.69 19.55 35.63 1.00 2189 2013
## y0[17] 15.93 15.79 5.68 5.83 6.78 25.35 1.00 1751 1656
## y0[18] 7.99 7.92 5.07 4.86 -0.37 16.54 1.00 2014 1777
## y0[19] 8.14 7.93 4.95 4.78 0.04 16.28 1.00 2163 1917
## y0[20] 20.05 20.02 5.59 5.65 10.75 29.32 1.00 2254 2014
## y0[21] 5.44 5.45 6.50 6.52 -5.18 16.30 1.00 1376 1461
## y0[22] 20.85 20.90 5.21 5.19 12.21 29.18 1.00 2065 1931
## y0[23] 24.40 24.47 5.35 5.17 15.66 33.13 1.00 2176 1973
## y0[24] 15.23 15.29 6.44 6.61 5.02 25.62 1.00 2036 1729
## y0[25] 35.58 35.47 4.79 4.71 27.63 43.48 1.00 2053 1715
## y0[26] 16.57 16.49 4.68 4.51 9.08 24.25 1.00 2123 1931
## y0[27] 47.23 47.24 4.80 4.67 39.50 55.17 1.00 2017 1813
## y0[28] 22.33 22.35 4.66 4.61 14.45 29.91 1.00 1932 1880
## y0[29] 37.46 37.37 6.23 6.10 27.43 47.88 1.00 2015 2046
## y0[30] 38.62 38.75 4.79 4.44 30.58 46.50 1.00 1920 1845
## y0[31] 18.45 18.44 5.51 5.22 9.52 27.59 1.00 1689 1728
## y0[32] 35.11 35.01 4.93 4.73 26.92 43.24 1.00 1873 1894
## y0[33] 5.62 5.80 5.20 5.06 -2.95 14.01 1.00 2202 1836
## y0[34] 53.39 53.32 5.22 5.10 44.92 62.26 1.00 2037 2056
## y0[35] 31.20 31.30 4.91 4.90 23.21 39.19 1.00 1814 1932
## y0[36] 6.11 6.09 6.42 6.24 -3.85 16.83 1.00 1800 1731
## y0[37] 22.54 22.51 5.18 4.91 13.99 30.85 1.00 1950 1779
## y0[38] 39.06 38.97 4.76 4.62 31.28 46.82 1.00 1975 2010
## y0[39] 33.49 33.51 4.76 4.86 25.52 41.37 1.00 2051 1917
## y0[40] 31.93 32.04 5.33 5.23 23.14 40.63 1.00 2008 1814
## y0[41] 32.52 32.43 4.72 4.69 25.07 40.61 1.00 2106 1865
## y0[42] 51.44 51.36 4.82 4.84 43.64 59.17 1.00 1932 1805
## y0[43] 11.54 11.51 4.79 4.69 3.71 19.50 1.00 2166 1568
## y0[44] 27.43 27.42 4.80 4.63 19.66 35.34 1.00 1900 1821
## y0[45] 33.77 33.85 4.98 4.80 25.50 41.94 1.00 2077 1946
## y0[46] 37.66 37.64 5.77 5.79 28.64 46.88 1.00 1867 1618
## y0[47] 37.23 37.29 4.69 4.62 29.58 44.94 1.00 1820 2092
## y0[48] 33.06 32.94 4.94 4.70 25.15 41.17 1.00 2165 1787
## y0[49] 7.17 7.23 5.05 4.96 -1.29 15.50 1.00 2082 2057
## y0[50] 43.94 44.00 5.22 5.11 35.67 52.61 1.00 2167 1865
m0=mcmc$draws('m0')
smm0=summary(m0)
y0=mcmc$draws('y0')
smy0=summary(y0)
grid.arrange(
qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))
par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')
xy=tibble(x=x,c=c,m=smm0$median,ml5=smm0$q5,mu5=smm0$q95,yl5=smy0$q5,yu5=smy0$q95)
qplot(x,y,shape=c,size=I(2),col=I('red'))+
scale_shape_manual(values=1:k) +
#geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
#geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
#geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
#geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=m),xy,col='black')
qplot(x,y,facets=~c,size=I(2),col=I('red'))+
geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=m),xy,col='black')
estimate as class have relation each other
all b0l,b1l have a distribution
class c=1~k
b0c~N(b00,sb0)
b1c~N(b10,sb1)
yc~N(b0c+b1c*x,s)ex12-2.stan
data{
int N;
int K;
vector[N] x;
vector[N] y;
array[N] int c;
}
parameters{
real b00;
real<lower=0,upper=100> sb0;
vector[K] b0;
real b10;
real<lower=0,upper=100> sb1;
vector[K] b1;
real<lower=0,upper=100> s;
}
model{
b0~normal(b00,sb0);
b1~normal(b10,sb1);
for(i in 1:N)
y[i]~normal(b0[c[i]]+b1[c[i]]*x[i],s);
}
generated quantities{
vector[N] m0;
vector[N] y0;
for(i in 1:N){
m0[i]=b0[c[i]]+b1[c[i]]*x[i];
y0[i]=normal_rng(m0[i],s);
}
}mdl=cmdstan_model('./ex12-2.stan')
data=list(N=n,K=k,x=x,y=y,c=factor(c))
mle=mdl$optimize(data=data)## Initial log joint probability = -272.182
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 99 -20.0657 0.276822 604236 0.5672 0.5672 141
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 121 -16.6215 8.53733e-05 290.799 1 1 167
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.try(print(mle))## variable estimate
## lp__ -16.62
## b00 1.02
## sb0 0.00
## b0[1] 1.02
## b0[2] 1.02
## b0[3] 1.02
## b0[4] 1.02
## b0[5] 1.02
## b0[6] 1.02
## b0[7] 1.02
## b0[8] 1.02
## b0[9] 1.02
## b10 2.30
## sb1 0.20
## b1[1] 2.72
## b1[2] 2.05
## b1[3] 2.96
## b1[4] 2.03
## b1[5] 2.36
## b1[6] 2.61
## b1[7] 1.52
## b1[8] 1.37
## b1[9] 1.93
## s 13.83
## m0[1] 11.32
## m0[2] 11.26
## m0[3] 35.45
## m0[4] 15.84
## m0[5] 47.05
## m0[6] 41.17
## m0[7] 16.74
## m0[8] 11.49
## m0[9] 49.95
## m0[10] 33.40
## m0[11] 2.87
## m0[12] 51.54
## m0[13] 46.88
## m0[14] 20.19
## m0[15] 24.58
## m0[16] 25.20
## m0[17] 11.60
## m0[18] 10.53
## m0[19] 6.35
## m0[20] 9.88
## m0[21] 5.74
## m0[22] 21.82
## m0[23] 25.08
## m0[24] 7.50
## m0[25] 25.95
## m0[26] 16.04
## m0[27] 46.18
## m0[28] 22.79
## m0[29] 24.57
## m0[30] 37.63
## m0[31] 13.13
## m0[32] 25.32
## m0[33] 8.24
## m0[34] 56.87
## m0[35] 31.72
## m0[36] 13.93
## m0[37] 14.49
## m0[38] 38.11
## m0[39] 32.08
## m0[40] 23.04
## m0[41] 30.95
## m0[42] 53.55
## m0[43] 10.35
## m0[44] 28.51
## m0[45] 36.75
## m0[46] 28.25
## m0[47] 28.75
## m0[48] 21.51
## m0[49] 9.66
## m0[50] 47.50
## y0[1] -6.47
## y0[2] -8.59
## y0[3] 38.94
## y0[4] 24.68
## y0[5] 35.80
## y0[6] 45.48
## y0[7] 22.54
## y0[8] 9.63
## y0[9] 51.73
## y0[10] 44.26
## y0[11] -22.82
## y0[12] 54.11
## y0[13] 60.51
## y0[14] 22.77
## y0[15] 33.30
## y0[16] 23.19
## y0[17] 19.22
## y0[18] 18.87
## y0[19] 16.99
## y0[20] 21.14
## y0[21] -4.88
## y0[22] 32.99
## y0[23] 25.37
## y0[24] 16.04
## y0[25] 23.39
## y0[26] 11.39
## y0[27] 45.55
## y0[28] 16.25
## y0[29] 30.36
## y0[30] 41.61
## y0[31] 21.84
## y0[32] 27.37
## y0[33] 8.73
## y0[34] 40.02
## y0[35] 51.81
## y0[36] 38.49
## y0[37] 41.32
## y0[38] 40.55
## y0[39] 26.07
## y0[40] 49.59
## y0[41] 47.88
## y0[42] 47.04
## y0[43] 2.00
## y0[44] 19.31
## y0[45] 28.59
## y0[46] 7.95
## y0[47] 45.51
## y0[48] 28.33
## y0[49] 27.21
## y0[50] 39.47mcmc=goMCMC(mdl,data,wrm=1000,smp=2000)## Running MCMC with 4 parallel chains...
##
## Chain 1 Iteration: 1 / 3000 [ 0%] (Warmup)
## Chain 2 Iteration: 1 / 3000 [ 0%] (Warmup)
## Chain 3 Iteration: 1 / 3000 [ 0%] (Warmup)
## Chain 4 Iteration: 1 / 3000 [ 0%] (Warmup)
## Chain 1 Iteration: 1000 / 3000 [ 33%] (Warmup)
## Chain 1 Iteration: 1001 / 3000 [ 33%] (Sampling)
## Chain 2 Iteration: 1000 / 3000 [ 33%] (Warmup)
## Chain 2 Iteration: 1001 / 3000 [ 33%] (Sampling)
## Chain 3 Iteration: 1000 / 3000 [ 33%] (Warmup)
## Chain 3 Iteration: 1001 / 3000 [ 33%] (Sampling)
## Chain 4 Iteration: 1000 / 3000 [ 33%] (Warmup)
## Chain 4 Iteration: 1001 / 3000 [ 33%] (Sampling)
## Chain 1 Iteration: 2000 / 3000 [ 66%] (Sampling)
## Chain 2 Iteration: 2000 / 3000 [ 66%] (Sampling)
## Chain 3 Iteration: 2000 / 3000 [ 66%] (Sampling)
## Chain 4 Iteration: 2000 / 3000 [ 66%] (Sampling)
## Chain 1 Iteration: 3000 / 3000 [100%] (Sampling)
## Chain 1 finished in 1.2 seconds.
## Chain 2 Iteration: 3000 / 3000 [100%] (Sampling)
## Chain 3 Iteration: 3000 / 3000 [100%] (Sampling)
## Chain 4 Iteration: 3000 / 3000 [100%] (Sampling)
## Chain 2 finished in 1.3 seconds.
## Chain 3 finished in 1.3 seconds.
## Chain 4 finished in 1.3 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 1.3 seconds.
## Total execution time: 1.4 seconds.seeMCMC(mcmc,'m')## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -116.83 -116.60 5.04 4.69 -125.45 -109.01 1.00 830 355
## b00 7.89 7.90 3.80 3.40 1.85 13.90 1.00 4118 2995
## sb0 9.57 8.88 3.64 3.04 5.10 16.51 1.00 4730 4995
## b0[1] 16.28 16.44 7.04 6.97 4.60 27.51 1.00 3612 1665
## b0[2] 10.98 10.84 4.12 4.08 4.42 17.88 1.00 2449 2785
## b0[3] 19.46 19.57 2.99 2.94 14.49 24.23 1.00 6108 4781
## b0[4] 0.36 0.24 3.87 3.81 -5.77 6.80 1.00 5684 3885
## b0[5] 7.24 7.28 3.95 3.96 0.78 13.71 1.00 2590 1211
## b0[6] 4.81 4.78 2.31 2.28 1.07 8.81 1.00 2195 672
## b0[7] 7.28 7.27 4.34 4.18 0.22 14.45 1.00 6415 5961
## b0[8] -2.17 -2.27 3.71 3.70 -8.04 4.04 1.00 3033 4060
## b0[9] 6.77 6.98 4.48 4.34 -1.02 13.66 1.00 4506 4735
## b10 1.82 1.82 0.26 0.22 1.42 2.23 1.00 5230 5290
## sb1 0.53 0.47 0.32 0.26 0.13 1.10 1.00 749 295
## b1[1] 2.17 2.10 0.55 0.49 1.40 3.18 1.00 4375 4999
## b1[2] 1.37 1.38 0.29 0.30 0.89 1.84 1.00 1515 715
## b1[3] 1.83 1.83 0.24 0.23 1.45 2.22 1.00 6103 4864
## b1[4] 1.52 1.61 0.59 0.49 0.39 2.33 1.00 5070 4319
## b1[5] 2.01 2.01 0.27 0.27 1.58 2.45 1.00 4253 5535
## b1[6] 2.13 2.13 0.23 0.23 1.75 2.50 1.00 2406 1383
## b1[7] 1.72 1.73 0.31 0.30 1.20 2.22 1.00 6744 5543
## b1[8] 1.54 1.56 0.31 0.31 1.00 2.01 1.00 3450 5350
## b1[9] 2.16 2.10 0.49 0.45 1.47 3.04 1.00 4213 4278
## s 4.49 4.45 0.57 0.55 3.65 5.49 1.00 4858 5462
## m0[1] 13.22 13.20 1.69 1.65 10.50 16.14 1.00 2866 951
## m0[2] 15.96 15.94 2.93 2.95 11.21 20.76 1.00 2379 1251
## m0[3] 34.06 34.00 2.21 2.18 30.50 37.85 1.00 2083 926
## m0[4] 28.64 28.71 2.04 1.98 25.24 31.85 1.00 6095 6002
## m0[5] 46.43 46.41 2.27 2.22 42.78 50.17 1.00 6814 5594
## m0[6] 37.90 37.87 2.75 2.74 33.43 42.63 1.00 1694 776
## m0[7] 15.47 15.45 1.71 1.68 12.66 18.27 1.00 7277 4665
## m0[8] 25.94 26.04 2.30 2.25 22.14 29.58 1.00 6106 5587
## m0[9] 49.79 49.78 1.95 1.95 46.62 52.96 1.00 6203 5561
## m0[10] 42.18 42.30 4.36 4.27 34.87 49.18 1.00 5719 4167
## m0[11] 6.31 6.28 2.18 2.16 2.77 10.10 1.00 2260 731
## m0[12] 46.06 46.08 2.92 2.89 41.27 50.86 1.00 4855 4103
## m0[13] 46.28 46.26 2.25 2.20 42.66 50.00 1.00 6816 5533
## m0[14] 19.34 19.33 1.98 1.97 16.07 22.56 1.00 6974 4742
## m0[15] 24.04 24.02 1.42 1.40 21.73 26.35 1.00 7464 5450
## m0[16] 27.82 27.82 1.83 1.83 24.88 30.83 1.00 2318 1381
## m0[17] 18.61 18.62 2.74 2.74 14.07 22.99 1.00 5611 3492
## m0[18] 8.50 8.49 2.05 2.02 5.15 11.88 1.00 4283 4903
## m0[19] 9.16 9.14 1.96 1.93 6.00 12.53 1.00 2424 890
## m0[20] 16.91 16.82 3.03 3.01 12.03 21.98 1.00 3155 4164
## m0[21] 3.88 3.84 3.22 3.20 -1.39 9.30 1.00 6481 4079
## m0[22] 21.17 21.19 2.20 2.20 17.52 24.70 1.00 6348 4749
## m0[23] 24.83 24.90 2.72 2.73 20.24 29.19 1.00 5634 5314
## m0[24] 14.63 14.64 3.36 3.27 9.18 20.15 1.00 6068 6436
## m0[25] 34.91 34.93 1.58 1.56 32.28 37.45 1.00 6085 5966
## m0[26] 17.07 17.06 1.50 1.48 14.65 19.68 1.00 3864 1208
## m0[27] 47.45 47.45 1.76 1.76 44.56 50.34 1.00 6186 4117
## m0[28] 22.58 22.59 1.41 1.38 20.29 24.88 1.00 6814 6080
## m0[29] 33.12 33.04 3.42 3.33 27.57 38.89 1.00 5424 3682
## m0[30] 38.41 38.40 1.70 1.67 35.68 41.19 1.00 5188 4868
## m0[31] 20.32 20.36 2.62 2.62 16.00 24.54 1.00 6011 3603
## m0[32] 34.52 34.54 1.60 1.57 31.86 37.08 1.00 6094 6039
## m0[33] 5.93 5.92 2.38 2.37 2.04 9.90 1.00 3740 5018
## m0[34] 54.08 54.07 2.35 2.33 50.28 57.91 1.00 6288 5774
## m0[35] 31.56 31.54 1.96 1.93 28.36 34.80 1.00 2776 1263
## m0[36] 10.02 10.01 3.38 3.30 4.45 15.64 1.00 6434 5095
## m0[37] 20.00 19.99 2.53 2.50 15.85 24.22 1.00 4018 5127
## m0[38] 38.82 38.80 1.72 1.68 36.06 41.64 1.00 5346 4916
## m0[39] 33.68 33.65 1.62 1.61 31.09 36.38 1.00 3026 1346
## m0[40] 32.27 32.28 2.75 2.68 27.80 36.72 1.00 6690 4214
## m0[41] 32.72 32.69 1.63 1.61 30.10 35.44 1.00 2837 1350
## m0[42] 52.02 51.99 2.15 2.15 48.56 55.52 1.00 6250 5641
## m0[43] 12.43 12.42 1.74 1.70 9.63 15.45 1.00 2748 953
## m0[44] 27.25 27.26 1.51 1.48 24.78 29.72 1.00 8183 5133
## m0[45] 34.93 34.88 2.32 2.28 31.20 38.88 1.00 1949 893
## m0[46] 38.17 38.22 3.28 3.21 32.79 43.49 1.00 6706 4179
## m0[47] 36.64 36.68 1.51 1.49 34.13 39.05 1.00 6145 6094
## m0[48] 32.16 32.20 1.75 1.71 29.21 34.93 1.00 6068 6292
## m0[49] 7.52 7.51 2.17 2.14 3.96 11.10 1.00 4055 4767
## m0[50] 42.76 42.77 2.62 2.60 38.47 47.04 1.00 5359 4143
## y0[1] 13.28 13.20 4.81 4.67 5.53 21.28 1.00 8258 8144
## y0[2] 15.90 15.90 5.35 5.22 6.97 24.77 1.00 6654 7256
## y0[3] 34.05 34.06 4.99 4.98 25.97 42.21 1.00 7190 7706
## y0[4] 28.60 28.58 4.93 4.79 20.68 36.75 1.00 7469 7821
## y0[5] 46.49 46.52 5.05 4.99 38.15 54.72 1.00 7663 7126
## y0[6] 37.96 37.91 5.22 5.16 29.51 46.68 1.00 6884 6871
## y0[7] 15.50 15.44 4.86 4.80 7.56 23.44 1.00 8019 7363
## y0[8] 26.01 26.02 5.00 4.91 17.86 34.15 1.00 7951 7672
## y0[9] 49.84 49.85 4.89 4.80 41.82 57.77 1.00 7868 7973
## y0[10] 42.21 42.17 6.29 6.19 32.03 52.67 1.00 8362 7510
## y0[11] 6.31 6.23 5.01 4.87 -1.83 14.60 1.00 7173 7346
## y0[12] 46.16 46.16 5.34 5.29 37.42 55.04 1.00 6264 7520
## y0[13] 46.34 46.39 5.09 5.06 38.05 54.77 1.00 7682 7210
## y0[14] 19.29 19.27 4.86 4.77 11.36 27.25 1.00 6412 7675
## y0[15] 23.93 23.94 4.79 4.67 16.07 31.80 1.00 8057 7823
## y0[16] 27.76 27.79 4.95 4.90 19.70 35.96 1.00 6546 6999
## y0[17] 18.60 18.59 5.35 5.22 9.82 27.38 1.00 7258 7555
## y0[18] 8.45 8.48 4.94 4.86 0.17 16.63 1.00 7126 7178
## y0[19] 9.12 9.12 4.95 4.91 0.94 17.26 1.00 7264 7303
## y0[20] 16.95 16.99 5.46 5.35 7.92 25.85 1.00 5686 6882
## y0[21] 3.88 3.83 5.59 5.47 -5.29 12.99 1.00 7259 7484
## y0[22] 21.20 21.27 4.95 4.88 13.11 29.30 1.00 7226 7259
## y0[23] 24.78 24.81 5.27 5.14 16.08 33.45 1.00 6567 7769
## y0[24] 14.64 14.58 5.65 5.60 5.52 24.22 1.00 7575 7449
## y0[25] 34.94 34.95 4.90 4.83 26.84 42.95 1.00 7661 7374
## y0[26] 17.12 17.10 4.87 4.72 9.14 24.88 1.00 7516 7726
## y0[27] 47.50 47.47 4.83 4.74 39.51 55.40 1.00 7561 7589
## y0[28] 22.57 22.63 4.75 4.69 14.79 30.22 1.00 7934 7486
## y0[29] 33.06 33.06 5.73 5.65 23.71 42.46 1.00 7563 7209
## y0[30] 38.42 38.39 4.83 4.73 30.45 46.38 1.00 7657 7721
## y0[31] 20.35 20.32 5.26 5.14 11.78 29.00 1.00 8028 7635
## y0[32] 34.55 34.57 4.79 4.61 26.64 42.42 1.00 7674 7872
## y0[33] 5.97 5.92 5.22 5.11 -2.65 14.57 1.00 7727 7931
## y0[34] 54.18 54.14 5.10 4.92 45.70 62.58 1.00 7647 7877
## y0[35] 31.56 31.57 4.89 4.69 23.46 39.64 1.00 6710 7857
## y0[36] 10.06 10.03 5.67 5.69 0.69 19.38 1.00 8136 6913
## y0[37] 19.99 19.99 5.18 5.15 11.39 28.35 1.00 7407 7749
## y0[38] 38.80 38.77 4.78 4.61 30.87 46.63 1.00 7407 7516
## y0[39] 33.67 33.68 4.80 4.69 25.75 41.66 1.00 7398 7442
## y0[40] 32.19 32.24 5.28 5.23 23.56 40.85 1.00 8265 7337
## y0[41] 32.71 32.66 4.81 4.75 24.95 40.65 1.00 7946 7759
## y0[42] 51.97 51.95 4.94 4.84 43.82 60.08 1.00 8211 8027
## y0[43] 12.47 12.45 4.82 4.71 4.64 20.34 1.00 7399 7045
## y0[44] 27.41 27.45 4.76 4.69 19.69 35.05 1.00 8157 7715
## y0[45] 34.97 35.00 5.06 4.96 26.68 43.26 1.00 6921 7626
## y0[46] 38.20 38.23 5.65 5.55 28.85 47.50 1.00 7330 6932
## y0[47] 36.74 36.71 4.73 4.56 28.84 44.56 1.00 7638 7483
## y0[48] 32.11 32.14 4.85 4.70 24.17 40.09 1.00 7992 7425
## y0[49] 7.51 7.50 5.05 4.86 -0.82 15.67 1.00 7264 7583
## y0[50] 42.72 42.71 5.20 5.13 34.21 51.16 1.00 7826 7712
m0=mcmc$draws('m0')
smm0=summary(m0)
y0=mcmc$draws('y0')
smy0=summary(y0)
grid.arrange(
qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))
par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')
xy=tibble(x=x,c=c,m=smm0$median,ml5=smm0$q5,mu5=smm0$q95,yl5=smy0$q5,yu5=smy0$q95)
qplot(x,y,shape=c,size=I(2),col=I('red'))+
scale_shape_manual(values=1:k) +
#geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
#geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
#geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
#geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=m),xy,col='black')
qplot(x,y,facets=~c,size=I(2),col=I('red'))+
geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=m),xy,col='black')
ex13
generalized linear mixed model
(X,y)i=1-n
b[b0,b1,...]
linear model y~N(Xb,s)
generalized linear model y~dist.(m=link(Xb),s)
fixed effect b0, b1,...
individual random effect b0+r0i r0~N(0,sr0), b1+r1i r1~N(0,sr1),...
class c=1-k
class effect b0+r0c r0~N(0,sr0), b1+r1c r1~N(0,sr1),...
for y=dist.(m,s)
random intercept model m=b0+r0i+b1*x, m=b0+r0c+b1*x
random coefficient model m=b0+(b1+r1i)*x, m=b0+(b1+r1c)*x
mixed model m=b0+r0i+(b1+r1i)*x, m=b0+r0c+(b1+r1c)*x
note
@ yi=b0+b1*xi+r0i is not useful, ri is included in s
@ yi=b0+b1*xi+r0c is useful, repeated observation can be treated by class effect
@ when appling Poisson dist.(y is larger than 0, error is larger at large x),
but variance is larger than mean (over dispersion),
random intercept model is usefuln=20
x=runif(n,0,10)
r0=rnorm(n,0,1)
y=rpois(n,exp(2+0.1*x+r0))
qplot(x,y)
ex13-0.stan
data{
int N;
vector[N] x;
array[N] int y;
}
parameters{
real b0;
real b1;
}
model{
y~poisson_log(b0+b1*x);
}
generated quantities{
vector[N] m0;
array[N] int y0;
for(i in 1:N){
m0[i]=b0+b1*x[i];
y0[i]=poisson_log_rng(m0[i]);
}
}mdl=cmdstan_model('./ex13-0.stan')
data=list(N=n,x=x,y=y)
mle=mdl$optimize(data=data)## Initial log joint probability = 274.414
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 11 318.903 0.000798363 0.0165867 0.9917 0.9917 14
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ 318.90
## b0 1.90
## b1 0.10
## m0[1] 2.68
## m0[2] 1.91
## m0[3] 2.31
## m0[4] 2.12
## m0[5] 2.71
## m0[6] 2.05
## m0[7] 2.61
## m0[8] 2.81
## m0[9] 2.81
## m0[10] 2.12
## m0[11] 2.18
## m0[12] 2.50
## m0[13] 2.45
## m0[14] 2.05
## m0[15] 2.04
## m0[16] 2.77
## m0[17] 2.33
## m0[18] 2.37
## m0[19] 2.29
## m0[20] 2.14
## y0[1] 13.00
## y0[2] 8.00
## y0[3] 17.00
## y0[4] 5.00
## y0[5] 9.00
## y0[6] 3.00
## y0[7] 12.00
## y0[8] 18.00
## y0[9] 22.00
## y0[10] 12.00
## y0[11] 9.00
## y0[12] 13.00
## y0[13] 15.00
## y0[14] 12.00
## y0[15] 9.00
## y0[16] 14.00
## y0[17] 9.00
## y0[18] 9.00
## y0[19] 9.00
## y0[20] 11.00mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,'m')## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ 318.01 318.31 0.94 0.63 316.16 318.86 1.00 531 603
## b0 1.90 1.89 0.13 0.13 1.69 2.11 1.01 458 367
## b1 0.10 0.10 0.02 0.02 0.06 0.13 1.01 568 515
## m0[1] 2.68 2.68 0.09 0.08 2.53 2.82 1.00 1678 1264
## m0[2] 1.91 1.91 0.13 0.13 1.70 2.12 1.01 458 372
## m0[3] 2.31 2.31 0.07 0.07 2.19 2.42 1.01 631 831
## m0[4] 2.12 2.12 0.09 0.09 1.97 2.28 1.01 473 387
## m0[5] 2.71 2.71 0.09 0.09 2.56 2.85 1.00 1561 1264
## m0[6] 2.06 2.05 0.10 0.10 1.89 2.23 1.01 461 375
## m0[7] 2.60 2.60 0.08 0.07 2.48 2.72 1.00 1970 1301
## m0[8] 2.80 2.81 0.11 0.11 2.63 2.98 1.00 1254 1313
## m0[9] 2.80 2.80 0.11 0.11 2.62 2.97 1.00 1269 1296
## m0[10] 2.12 2.12 0.09 0.09 1.97 2.28 1.01 473 384
## m0[11] 2.18 2.19 0.08 0.08 2.05 2.33 1.01 497 499
## m0[12] 2.49 2.50 0.07 0.06 2.38 2.60 1.00 1556 1384
## m0[13] 2.44 2.44 0.06 0.06 2.34 2.55 1.01 1163 1233
## m0[14] 2.05 2.05 0.11 0.11 1.88 2.22 1.01 460 375
## m0[15] 2.04 2.04 0.11 0.11 1.87 2.22 1.01 460 375
## m0[16] 2.76 2.76 0.10 0.10 2.60 2.92 1.00 1363 1367
## m0[17] 2.33 2.33 0.07 0.06 2.22 2.44 1.01 677 929
## m0[18] 2.36 2.36 0.07 0.06 2.25 2.47 1.01 778 978
## m0[19] 2.29 2.29 0.07 0.07 2.18 2.41 1.01 605 804
## m0[20] 2.14 2.14 0.09 0.09 1.99 2.29 1.01 478 426
## y0[1] 14.55 14.00 3.99 4.45 8.00 21.00 1.00 2075 1926
## y0[2] 6.80 7.00 2.74 2.97 3.00 12.00 1.00 1637 1760
## y0[3] 10.01 10.00 3.18 2.97 5.00 15.00 1.00 1803 1814
## y0[4] 8.32 8.00 3.02 2.97 4.00 13.05 1.00 1583 1876
## y0[5] 14.92 15.00 4.06 4.45 9.00 22.00 1.00 1872 1915
## y0[6] 7.87 8.00 2.94 2.97 3.00 13.00 1.00 1536 1868
## y0[7] 13.41 13.00 3.80 2.97 8.00 20.00 1.00 2108 1852
## y0[8] 16.56 16.00 4.42 4.45 10.00 24.00 1.00 2048 1983
## y0[9] 16.45 16.00 4.45 4.45 10.00 24.00 1.00 1930 1988
## y0[10] 8.41 8.00 2.99 2.97 4.00 14.00 1.00 1885 1877
## y0[11] 9.06 9.00 3.08 2.97 4.00 14.00 1.00 1721 1751
## y0[12] 12.16 12.00 3.68 2.97 6.00 18.00 1.00 2089 1828
## y0[13] 11.46 11.00 3.39 2.97 6.00 17.00 1.00 2073 2037
## y0[14] 7.78 8.00 2.87 2.97 3.00 13.00 1.00 1746 2097
## y0[15] 7.74 8.00 2.89 2.97 3.00 13.00 1.00 1683 1964
## y0[16] 16.04 16.00 4.36 4.45 9.00 24.00 1.00 1944 1867
## y0[17] 10.21 10.00 3.35 2.97 5.00 16.00 1.00 1804 1870
## y0[18] 10.79 11.00 3.44 2.97 5.00 17.00 1.00 1887 1949
## y0[19] 9.87 10.00 3.10 2.97 5.00 15.00 1.00 2014 1997
## y0[20] 8.51 8.00 3.04 2.97 4.00 14.00 1.00 1736 1890
m0=mcmc$draws('m0')
smm0=summary(m0)
y0=mcmc$draws('y0')
smy0=summary(y0)
grid.arrange(
qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))
par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')
xy=tibble(x=x,m=smm0$median,ml5=smm0$q5,mu5=smm0$q95,yl5=smy0$q5,yu5=smy0$q95)
qplot(x,y,size=I(2),col=I('red'))+
scale_shape_manual(values=1:k) +
geom_line(aes(x=x,y=exp(ml5)),xy,col='darkgray')+
geom_line(aes(x=x,y=exp(mu5)),xy,col='darkgray')+
geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=exp(m)),xy,col='black')
ex13-1.stan
data{
int N;
vector[N] x;
array[N] int y;
}
parameters{
real b0;
real b1;
real<lower=0> sr0;
vector[N] r0;
}
model{
r0~normal(0,sr0);
for(i in 1:N)
y~poisson_log(b0+r0[i]+b1*x);
}
generated quantities{
vector[N] m0;
array[N] int y0;
for(i in 1:N){
m0[i]=b0+r0[i]+b1*x[i];
y0[i]=poisson_log_rng(m0[i]);
}
}mdl=cmdstan_model('./ex13-1.stan')
data=list(N=n,x=x,y=y)
mle=mdl$optimize(data=data)## Initial log joint probability = -3847.38
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 99 6539.77 0.0415256 8906.86 0.2214 0.2214 146
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 152 6658.97 0.00244894 1.88547e+06 0.8696 0.8696 246
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ 6658.97
## b0 1.85
## b1 0.14
## sr0 0.00
## r0[1] 0.00
## r0[2] 0.00
## r0[3] 0.00
## r0[4] 0.00
## r0[5] 0.00
## r0[6] 0.00
## r0[7] 0.00
## r0[8] 0.00
## r0[9] 0.00
## r0[10] 0.00
## r0[11] 0.00
## r0[12] 0.00
## r0[13] 0.00
## r0[14] 0.00
## r0[15] 0.00
## r0[16] 0.00
## r0[17] 0.00
## r0[18] 0.00
## r0[19] 0.00
## r0[20] 0.00
## m0[1] 2.98
## m0[2] 1.87
## m0[3] 2.44
## m0[4] 2.17
## m0[5] 3.03
## m0[6] 2.08
## m0[7] 2.87
## m0[8] 3.17
## m0[9] 3.16
## m0[10] 2.18
## m0[11] 2.26
## m0[12] 2.72
## m0[13] 2.64
## m0[14] 2.07
## m0[15] 2.05
## m0[16] 3.11
## m0[17] 2.47
## m0[18] 2.53
## m0[19] 2.42
## m0[20] 2.20
## y0[1] 24.00
## y0[2] 4.00
## y0[3] 13.00
## y0[4] 10.00
## y0[5] 25.00
## y0[6] 5.00
## y0[7] 16.00
## y0[8] 22.00
## y0[9] 20.00
## y0[10] 12.00
## y0[11] 2.00
## y0[12] 14.00
## y0[13] 15.00
## y0[14] 6.00
## y0[15] 6.00
## y0[16] 20.00
## y0[17] 12.00
## y0[18] 14.00
## y0[19] 12.00
## y0[20] 15.00mcmc=goMCMC(mdl,data,wrm=500)## Running MCMC with 4 parallel chains...
##
## Chain 1 Iteration: 1 / 1000 [ 0%] (Warmup)
## Chain 2 Iteration: 1 / 1000 [ 0%] (Warmup)
## Chain 3 Iteration: 1 / 1000 [ 0%] (Warmup)
## Chain 4 Iteration: 1 / 1000 [ 0%] (Warmup)
## Chain 1 Iteration: 501 / 1000 [ 50%] (Sampling)
## Chain 2 Iteration: 501 / 1000 [ 50%] (Sampling)
## Chain 3 Iteration: 501 / 1000 [ 50%] (Sampling)
## Chain 4 Iteration: 501 / 1000 [ 50%] (Sampling)
## Chain 2 Iteration: 1000 / 1000 [100%] (Sampling)
## Chain 2 finished in 0.7 seconds.
## Chain 3 Iteration: 1000 / 1000 [100%] (Sampling)
## Chain 3 finished in 0.8 seconds.
## Chain 1 Iteration: 1000 / 1000 [100%] (Sampling)
## Chain 4 Iteration: 1000 / 1000 [100%] (Sampling)
## Chain 1 finished in 0.9 seconds.
## Chain 4 finished in 0.8 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.8 seconds.
## Total execution time: 1.0 seconds.seeMCMC(mcmc,'[m,r]')## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ 6454.08 6452.80 14.62 17.20 6431.54 6477.87 1.06 61 96
## b0 1.90 1.90 0.03 0.03 1.85 1.95 1.01 444 913
## b1 0.10 0.10 0.01 0.01 0.09 0.11 1.01 424 1034
## sr0 0.01 0.01 0.01 0.01 0.00 0.03 1.05 63 111
## r0[1] 0.00 0.00 0.01 0.01 -0.02 0.03 1.02 2753 857
## r0[2] 0.00 0.00 0.02 0.01 -0.03 0.02 1.01 3486 863
## r0[3] 0.00 0.00 0.02 0.01 -0.02 0.03 1.02 3180 802
## r0[4] 0.00 0.00 0.02 0.01 -0.02 0.02 1.02 4061 452
## r0[5] 0.00 0.00 0.01 0.01 -0.02 0.02 1.02 3127 838
## r0[6] 0.00 0.00 0.02 0.01 -0.03 0.02 1.02 2805 720
## r0[7] 0.00 0.00 0.02 0.01 -0.03 0.02 1.04 3817 799
## r0[8] 0.00 0.00 0.01 0.01 -0.02 0.02 1.03 3862 646
## r0[9] 0.00 0.00 0.02 0.01 -0.03 0.03 1.04 2789 860
## r0[10] 0.00 0.00 0.02 0.01 -0.03 0.03 1.02 4005 832
## r0[11] 0.00 0.00 0.02 0.01 -0.02 0.02 1.03 3282 733
## r0[12] 0.00 0.00 0.02 0.01 -0.03 0.03 1.02 3616 784
## r0[13] 0.00 0.00 0.02 0.01 -0.03 0.02 1.02 3440 630
## r0[14] 0.00 0.00 0.02 0.01 -0.03 0.03 1.03 3249 751
## r0[15] 0.00 0.00 0.02 0.01 -0.02 0.02 1.02 3095 919
## r0[16] 0.00 0.00 0.01 0.01 -0.02 0.02 1.01 2387 508
## r0[17] 0.00 0.00 0.02 0.01 -0.03 0.02 1.03 3014 750
## r0[18] 0.00 0.00 0.02 0.01 -0.02 0.02 1.02 3729 624
## r0[19] 0.00 0.00 0.01 0.01 -0.02 0.02 1.02 3096 934
## r0[20] 0.00 0.00 0.02 0.01 -0.02 0.02 1.03 4177 558
## m0[1] 2.68 2.68 0.02 0.02 2.64 2.72 1.00 1961 1407
## m0[2] 1.91 1.91 0.04 0.04 1.85 1.96 1.01 519 868
## m0[3] 2.31 2.31 0.02 0.02 2.27 2.34 1.00 1491 1371
## m0[4] 2.12 2.12 0.03 0.03 2.08 2.17 1.01 735 852
## m0[5] 2.71 2.71 0.02 0.02 2.67 2.75 1.00 1292 1573
## m0[6] 2.05 2.05 0.03 0.03 2.01 2.10 1.00 613 1058
## m0[7] 2.61 2.61 0.02 0.02 2.57 2.64 1.02 2856 1257
## m0[8] 2.81 2.81 0.03 0.03 2.77 2.86 1.00 1072 1543
## m0[9] 2.81 2.81 0.03 0.03 2.76 2.85 1.00 1062 1568
## m0[10] 2.12 2.12 0.03 0.03 2.08 2.17 1.00 811 1492
## m0[11] 2.18 2.18 0.03 0.03 2.14 2.22 1.00 930 1025
## m0[12] 2.50 2.50 0.02 0.02 2.46 2.53 1.01 2913 1250
## m0[13] 2.45 2.45 0.02 0.02 2.41 2.48 1.02 2630 1082
## m0[14] 2.05 2.05 0.03 0.03 2.00 2.10 1.01 734 1025
## m0[15] 2.04 2.04 0.03 0.03 1.99 2.09 1.00 555 1289
## m0[16] 2.77 2.77 0.03 0.03 2.73 2.81 1.00 1215 1528
## m0[17] 2.33 2.33 0.02 0.02 2.29 2.36 1.01 1604 1008
## m0[18] 2.37 2.37 0.02 0.02 2.33 2.40 1.02 1982 1217
## m0[19] 2.29 2.29 0.02 0.02 2.26 2.33 1.00 1130 1066
## m0[20] 2.14 2.14 0.03 0.03 2.09 2.18 1.00 816 1343
## y0[1] 14.71 15.00 3.81 4.45 9.00 21.00 1.00 1711 1469
## y0[2] 6.77 7.00 2.57 2.97 3.00 11.00 1.00 1734 1827
## y0[3] 9.99 10.00 3.16 2.97 5.00 16.00 1.00 2016 1872
## y0[4] 8.28 8.00 2.78 2.97 4.00 13.00 1.00 2116 1846
## y0[5] 15.10 15.00 3.95 4.45 9.00 22.00 1.00 1915 1944
## y0[6] 7.93 8.00 2.83 2.97 3.95 13.00 1.00 1944 1911
## y0[7] 13.46 13.00 3.58 2.97 8.00 20.00 1.00 2058 2014
## y0[8] 16.65 17.00 4.06 4.45 10.00 23.05 1.00 1976 2053
## y0[9] 16.54 16.00 4.09 4.45 10.00 23.00 1.00 1951 1958
## y0[10] 8.36 8.00 2.88 2.97 4.00 13.00 1.00 2056 1983
## y0[11] 8.83 9.00 2.97 2.97 4.00 14.00 1.00 2017 1821
## y0[12] 12.17 12.00 3.58 2.97 6.00 18.00 1.00 1853 2006
## y0[13] 11.50 11.00 3.45 2.97 6.00 18.00 1.00 1900 1983
## y0[14] 7.72 8.00 2.76 2.97 4.00 13.00 1.00 2114 1902
## y0[15] 7.82 8.00 2.84 2.97 3.00 13.00 1.00 2089 1977
## y0[16] 15.89 16.00 3.93 4.45 10.00 23.00 1.00 2060 1680
## y0[17] 10.32 10.00 3.24 2.97 5.00 16.00 1.00 1936 2068
## y0[18] 10.64 10.50 3.26 3.71 6.00 16.00 1.00 1922 1898
## y0[19] 9.99 10.00 3.19 2.97 5.00 15.00 1.00 1935 1899
## y0[20] 8.43 8.00 2.85 2.97 4.00 13.00 1.00 1866 1805
m0=mcmc$draws('m0')
smm0=summary(m0)
y0=mcmc$draws('y0')
smy0=summary(y0)
grid.arrange(
qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))
par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')
xy=tibble(x=x,m=smm0$median,ml5=smm0$q5,mu5=smm0$q95,yl5=smy0$q5,yu5=smy0$q95)
qplot(x,y,size=I(2),col=I('red'))+
scale_shape_manual(values=1:k) +
geom_line(aes(x=x,y=exp(ml5)),xy,col='darkgray')+
geom_line(aes(x=x,y=exp(mu5)),xy,col='darkgray')+
geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=exp(m)),xy,col='black')